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Abstract 

Period variability, quantified by the standard deviation (SD) of the cycle-to-cycle period, is 
investigated for noisy phase oscillators. We define the checkpoint phase as the beginning/end 
point of one oscillation cycle and derive an expression for the SD as a function of this phase. We 
find that the SD is dependent on the checkpoint phase only when oscillators are coupled. The 
applicability of our theory is verified using a realistic model. Our work clarifies the relationship 
between period variability and synchronization from which valuable information regarding coupling 
can be inferred. 
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Oscillators functioning as clocks, such as crystal oscillators spin torque oscillators 



-|5|, and circadian and heart pacemakers 



-|8|, play an important role in various systems. 



Although these clocks are subjected to various types of noise, including thermal, quantum, 
and molecular noise, they are required to perform temporally precise oscillations; i.e., os- 
cillations with only a small variability in the period (known as "period jitter" in electronic 
engineering 9|]). 

In many cases, it is sufficient for the clock to strike precisely at a specific time in each 
oscillation cycle, and thus a perfectly regular oscillation waveform is not needed. For cardiac 
pacemakers only the moment of stimulation is relevant. Experimental data regarding circa- 



dian activity in mice 



10| indicate that the variability in the period between each activity 



onset is smaller than that between each offset. Similar results have also been obtained in 



explant circadian pacemaker tissue (the suprachiasmatic nucleus, SON) [10|. These obser- 
vations suggest that the onset is more important than the offset in a circadian clock, which 
may be designed in such a way that the crucial moment is expressed with high precision. 

Remember that the definition of an oscillation period requires a fixed beginning/ end point 
for each oscillation cycle; hereafter referred to as the checkpoint (Fig. [1]). Although the 
average period does not depend on the particular choice of checkpoint, the period variability 
may be sensitive to the checkpoint. In order to clarify whether the checkpoint dependence 
in circadian activity is an artifact due to a technical problem in determining the onset and 
offset times or an essential property of the circadian clock, we need to investigate under 
what conditions the period variability is dependent on the checkpoint; this has received 
scant attention to date. 

Another important aspect of the period variability is its relationship to synchronization. 
A clock is commonly synchronized to its master clock such as in the case of the SCN in 
response to the daily variation of sunlight, and in peripheral clocks in response to the SCN. 
In addition, most biological clocks, including the SCN, cardiac pacemakers, and pacemakers 
in weakly electrical fish, are composed of a population of synchronized oscillators js, 7, 



111- 



It is known, both experimentally and theoretically, that period variability is reduced when 



the oscillators are coupled and synchronized 



Eicaiiy, tr 
6|,|3]l6 



16| . The question, therefore, arises as to 



whether the checkpoint dependence of the period variability is attributable to the interaction 
between oscillators. 

In this Letter, we discuss this checkpoint dependence for the case of coupled noisy phase 



oscillators. The period variability can be quantified using the standard deviation (SD) of the 
cycle-to-cycle period, and we show that although the SD is not dependent on the checkpoint 
in a single phase oscillator, it is dependent in a system of coupled phase oscillators; i.e., the 
checkpoint dependence results from the coupling effect. The SD is derived as a function of 
the checkpoint phase, which clarifies the relationship between the SD and synchronization. 
In particular, we find that in the case of diffusive coupling between oscillators, the checkpoint 
dependence of the SD has the same tendency as that of the synchronization: the SD is small 
when the oscillators are well synchronized. In other cases, however, the relationship is more 
complex. We also apply our theory to a realistic model of the electrical activity in a cell to 
demonstrate its validity. We believe that this is the first theoretical study to elucidate the 
existence of precise timing and its relationship with synchronization. 

To begin, we prove that the period variability is independent of the checkpoint in a 
single phase oscillator system. When a limit cycle oscillator is subjected to weak noise, its 



18| 



dynamics are well described by the following phase oscillator model 17 

f^=uj + z{e).fDm, (1) 

where 9 and u are the phase and natural frequency, respectively. The 27r-periodic function 
Z{9) is a phase sensitivity function, which quantifies the phase response of the oscillator to 
noise, and ^(t) denotes independent and identically distributed (i.i.d.) noise; each random 
variable ^{t) for all t obeys the same probability distribution and all are mutually indepen- 
dent. The positive constant D denotes the noise strength. Note that our proof below holds 
even if we permit u and the probability distribution of ^ to be 27r-periodic functions of 9: 
uj{9) and ^{t,9). 

The kth oscillation time of an oscillator, t^i^'^ , is defined as the time at which 9 passes 
through 27rk + 9cp (0 < 9cp < 2tt) for the first time [Fig. ^h)]. We define 9cp as the 
checkpoint phase. The kth. oscillation period At^''^ is defined as At^''^ = t^'"'" — t^!!;^, and the 
SD is defined as 



SD(e,p) = Vi^[(AC-r)2], (2) 

where -&[■■■] represents the statistical average over k, and r is the average period given by 
r = E[At^i^^^]. Note that E[- ■ ■] denotes both the statistical average taken over k and the 
ensemble average in the present paper, which are identical in the steady state. The system 
given by Eq. ([1]) is always in the steady state. 



To prove that the SD is independent of 9cp, we introduce two checkpoint phases denoted 
by a and f3 [Fig. 1(b)]. Since the processes a — j- /3 and /3 — j- a for any k are independent, 
we arrive at SD(a) = SD(/3) for any arbitrary checkpoint phases a and (3. A detailed proof 
is given in Appendix A. 




(b) 













FIG. 1. (color online), (a) An example of the time series of an oscillation. Periods are observed at 
two checkpoints, a and /?. (b) The corresponding checkpoint phases in the phase description. 

Next, we consider a pair of coupled phase oscillators subjected to noise. When limit cycle 
oscillators are weakl y co upled to each other and subjected to weak noise, the dynamics can 
be described by flH] 



(3) 



02 = 00 + kJ{92, 9i) + Z{92)VD^2{t), 

where 6i and k > are the phase of the oscillator i and the coupling strength, respectively. 
The i.i.d. noise ^j(t) satisfies E[^i{t)] = and E[^i{t)^j(t')] = 6ij6{t — t'). The 27r-periodic 
function J(x, y) describes the interaction between oscillators, which leads to synchronization. 
We assume that, in the absence of noise (-D = 0), the oscillators are synchronized in phase, 
i.e., Oi^2{t) -+ 4>{t) {t — oo), where 0(t) is a solution of 

(j){t) =uj + U{4>,(f)). (4) 

The necessary condition for the stability of in-phase synchrony for Z) = is provided below 
[see Eq. ( ITT]) ]. We also assume that u + Kj{(f),(f)) > for any for the coupled system to 
be oscillatory. 

Our particular interest is in the relationship between the SD [Eq. ([2])] and the synchro- 
nization of two oscillators. We thus introduce the following order parameter that measures 
the phase distance from the in-phase state: 



d{9, 



^[ll^l-^2|IV 
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(5) 



where E[x(t)]gj^=0^^_^ represents the average of Xk over k (where Xk is the vahie of x(t) taken 
when 6*1 passes through 27: k + 6cp for the first time), and ||6'i — ^2!! is the phase difference 
defined on the ring [— vr, vr). The phase distance d{6cp) is zero when the oscillators are 
completely synchronized in phase, and increases with the phase difference. 

As we demonstrate below, the relationship between SD(6'cp) and d{6cp) is qualitatively 
different for the two cases where J(0, 0) is (A) independent of and (B) dependent on 0. 
Cases (A) and (B) imply that given in Eq. (jl]) is independent of and dependent on 
0, respectively. Phase reduction theory indicates that it is appropriate to assume the form 
J(x, y) = z{x)G{x, y), where z{x) is the phase sensitivity function for the interaction G{x, y) 



17 



It is known that diffusive coupling between chemical oscillators and gap-junction 



coupling between cells yields J{x,y) = z{x){h{x) — h{y)), where h represents a chemical 
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concentration 
allows the form J 
Kuramoto model 



or membrane potential, which corresponds to case (A). Case (A) also 
[^y) = ~ y)i which has been employed in many models such as the 
however, we do not employ this form in the demonstration, since the 



term j{x—y) is derived as a result of averaging the interaction z{x)G{x, y) over one oscillation 



period 18|, and, by this approximation, the information about the 6cp dependence is lost. 



Many other types of coupling, such as J{x,y) = z{x)h{y) employed below, correspond to 



case (B) 



2l|. 



As an example of case (A), we consider z{6) = sin 6* for < 6* < tt, z{6) = for 
7T < 9 < 27r, and h(6) = cos6, and the following as an example of case (B): z{6) = — sin^ 



and h{e) = 1 + cos ^ (2l|] . We set Z{e) = l,u = 27t, VD = 0.03 x 27r, and ^i(O) = ^2(0) = 0, 
and assume ^i,2(^) to be white Gaussian noise. We integrate Eq. ([3]) using the Euler scheme 
with a time step of 5 x 10~^ for t = 0-10100 and discard the t = 0-100 data as transient. 

Using these examples, numerically obtained SD values for 61 are plotted as a function 
of 9cp in Fig. |2]^a) and (b). The results indicate clearly the existence of O^p dependence 
in both cases, which was absent in the single phase oscillator system. This dependence 
becomes stronger for larger k values. In contrast, for k ^ w, the dependence vanishes 



because J{x,y) is well approximated by j{x — y) [18|, and thus, the system effectively has 
rotational symmetry. The O^p value at which SD(^cp) assumes its minimum represents the 
most precise timing. 

The 6cp dependence of d{6 cp) for the two cases is shown in Fig. |2] (c) and (d). A 
comparison with SD(^cp) shows that the checkpoint phase maxima and minima of each 



K value coincide in the case of (A). Thus, the most precise timing is obtained when the 
oscillators are synchronized. By contrast, the ^cp dependence is considerably different in 
the case of (B). Therefore, we expect that nontrivial factors, apart from synchronization, 
influence the SD. We also examined several other functions, z{9), h{6), and Z{6), and found 
a similar relationship between SD(^cp) and d{6cp) (data not shown). 




FIG. 2. (color online). The ST){6cp)/T for case (A) and (B) is shown in (a) and (b), respectively, 
where the vertical scale is expressed as a percentage. The distance from in-phase synchronization, 
d{6cp), for case (A) and (B) is shown in (c) and (d), respectively. The points and lines are 
the numerical results of the simulation and analytical predictions given by Eqs. (jl6p and (|13p . 
respectively. 

We now derive an expression for the SD. The derivation consists of two steps: (i) calcu- 
lation of the phase diffusion a{9cp) [defined by Eq. (JTj)] with a linear approximation, and (ii) 
transformation from a{6cp) to SD(0cp)- Here, we employ the solution 0(t) of Eq. (jl]) with 
0(0) = and the time tcp is defined by 0(tcp) = ^cp- The oscillation period for D = is 
denoted by r; i.e., 0(tcp + t) = O^p + 27r. After a transient time, our system approaches the 
steady state, which is defined by the following equation for all \t': 

Pm - ^2||; ^1 = *) = Pm - ^2||; ^1 = * + 27r), (6) 

where -P(||^i — 02\\] di = ^) is the probability density function of the distance ||^i — ^2!! at 
61 = We assume that the system is in the steady state at t = 0. The ensemble we consider 
here is defined by the initial condition at t = tcp, 6*1 (tcp) = Ocp, and 6*2 (tcp) is distributed in 
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[^i(^cp) — TT, ^i(tcp) + tt) according to Eq. (Q. From this point, E[- ■ ■ ] represents the average 
taken over this ensemble. The phase diffusion a{6cp) is defined by 

a(eep)' = i?[(^i(tcp + r) - ^i(tep) - (7) 

We also assume that the noise intensity D is sufficiently small and that the other parameters 
and functions are of 0(1), so that the phase difference ||^i — 6*2!! is small in most cases in 
the steady state. 

To calculate the phase diffusion, we decompose 61^2 as Oi^2(t) = (pit) + Ai^2(^)- We then 
consider the time duration < t < 0(r), in which Ai 2(t) <^ 1 is expected in most cases 
because D '^1. Therefore, we can linearize Eq. ([3]). We define the two modes, X = A1 + A2 
and y = Ai — A2, which obey 

{X,Y) = Kfx,Ym)){X,Y)+^xAtA{t)), (8) 



where M<P) ^ fL.,., + f _ _ = fyi<P) - a.l..,., 8y 



dJ{<l>,<l>) f,jA\ = dJ\ _ dJ 



, and 

x=y=if> 

^x,y(t,0(t)) = VDZ{(f){t)){^i{t) ±6W)- Note that = for all in case (A). 

The solutions of Eq. ([8]) can be described as 

=exp [+KFx,Y{<l>m 
X |(X,F)(0) + ^ exp [-KFx,y(0(t'))]ex,y(t',0(O)^^^'} , (9) 

where Fx,y(0(t)) = /x,y(0(i'))'^^'- Furthermore, because fx{<P) = ^6 obtain 

^i. (II). (10, 

For = 0, we obtain Fy(27r) = (1/k) ln(y(r)/y(0)). Therefore, in the absence of noise, 
in-phase synchronization is stable if 

Fy(27r) = c < 0. (11) 

The correlations, E[X(t)^], E[Y(tY], and E[X(t)Y{t)] are given in Appendix B. Since 
Eq. (ED can be rewritten as P(| Ai- A2I; t) = P(| Ai- A2I; t+r), then E[Y{tf] = E[Y{t+Tf] 
holds approximately, leading to 

E[Ym = 2D "^P'^^'J r Zm')rexp[-2nFy{m)W. (12) 



In addition, because d{6cpY = E[{Ai(tcp) — A2(tcp))^] = E[Y(tcp)'^], we obtain 

d(^cp)' =exp[2KFy(^cp)] 

X lE[YiO)']+2D Z(0)2exp[-2«:Fy(0)]^rf0j , (13) 

which is generally 6'cp- dependent even if Z^cf)) is constant. 

Using these correlations and Eq. ( ITOl) . we obtain the following expression for the phase 
diffusion (See Appendix B) 

a{e,pf = E[(Ai(t,p + r) - Ai(tep))'] 

= Cimvf + C2d{e,pf, (14) 

where the Ci^2 are independent of ^cp and are given by Ci = y ^^^^ ^"^^ ^2 = 
(1 — exp[Kc])/2. The Ci term is an effective diffusion constant for the center of the two 
oscillators, which is half that of an uncoupled oscillator, and the C2 term is associated with 
the stability of the synchronization. 

To transform a{9cp) to SD(6'cp), we note that when the noise intensity is low, most of the 
trajectories of ^i(t) are very close to the unperturbed trajectory (pit) (see Appendix C). In 
such the following relation approximately holds true: 

rr(R \ 

'cp)- (15) 



SD(^, 



cp; 



The same approximation (but for constant 0) was employed in Ref. 16| and verified nu- 
merically. 

From Eqs. fIBlip and (1151) . we finally arrive at 



SD(a 



cp; 



Cl + Cj^^^. (16) 



cp; 



The analytical results given by Eqs. (fT6|) and (fT3|l are in excellent agreement with the nu- 
merical results (Fig. |2]). Although we have only discussed paired identical phase oscillators, 
our theory can easily be extended to other cases, e.g., globally coupled (all-to-all) identical 
oscillators or a periodically driven noisy oscillator. 

Equation (ITB|) shows that the periodicity of SD(^cp) is based on the synchronization 
d{9cp) and phase velocity 0(6'cp). For case (A), since <p{9cp) is constant, there is one-to-one 
correspondence between SD(6'cp) and d{6cp)] i.e., the most precise timing (6'™™) is the timing 
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at which the best synchronization is achieved. This was observed in Fig. |2](a) and (c), where 
6*™" = 7r/2 + 0{k~^) can be obtained from dd{6cp) /dOcp = 0. For case (B), however, the SD 
also depends on 0(6'cp); this is in contrast to that observed for the single phase oscillator 
system in which the phase velocity uj{6) does not contribute to the checkpoint dependence 
of the SD. Figures |21^b) and (d) showed that SD(^^cp) and d{6cp) are considerably different, 
which indicates the strong effect of in this particular example. Indeed, SD(6'cp) assumes 
its minimum around a maximum (piOcp) (^™" ~ Stt/S). 

To investigate whether Eq. (fT6ll holds for a more realistic model, we employ the FitzHugh- 
Nagumo model given by 

Vi = Vi{V, - a)(l - V,) -Wi+ + Kv{V2 - Vi), ^^^^ 
Wi = e{Vi - bWi) + KwiW2 - Wi), 
in which the second oscillator is described in a similar way. We fixed a = —0.1, b = 0.5, 
and e = 0.01. This system shows limit-cycle oscillations with a period of r ~ 126.5 when 
noise and coupling are absent. The white Gaussian noise ^i(t) has an intensity of 0.01. The 
interaction is diffusive, i.e., case (A), and we consider the following two types: V^-coupling 
{Ky = 0.01, Kw = 0) and ly-coupling {Ky = 0,Kw = 0.01). The phase 9 was defined 
properly (see Appendix D), and SD(6'cp) and d{6cp) were obtained numerically. Figure [3] 
shows that the 6'cp-dependence of the SD is different in the two cases, suggesting a significant 
effect from the coupling. We estimated the Ci and C2 values using Eq. (16) and the least- 
squares method under the condition that both cases have the same Ci value, resulting in 
Ci = 5.4, C^^^ = 0.20, and C^^^ = 0.48. In Fig. [3l we can see that the SD is described well 
by Eq. ( fT6l) using the fitted Ci and C2 values. This demonstrates that the theory is valid 
for this biological model. 

In many cases, only the SD measured at a functionally relevant checkpoint characterizes 
the performance of a clock. When designing a precise clock, we only have to reduce SD(6'cp) 
for a specific 9cp- Equation (fT6l) implies that SD(^cp) at a given ^cp decreases with decreasing 
d{Ocp) and increasing (j){6cp). Therefore, attractive coupling between oscillators should be 
activated around the functionally relevant timing point. In addition, in case (B), the phase 
velocity should be increased through coupling. 

Our theory enables us to infer the coupling timing or form by measuring SD(0cp) at 
several checkpoints. Although this is, in principle, possible with d{9cp), using SD(6'cp) has 
the added advantages that the SD can be measured from a single time series and that d{6cp) 
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FIG. 3. (color online). Validation of Eq. (116p in the FitzHugh-Nagumo model. Open symbols are 
the numerically obtained SD values. Filled symbols are the SD values evaluated from Eq. ()16p 
with the numerically obtained d values and fitting parameters C\ and C^- The triangles and circles 
represent the V- and VF-coupling cases. The plus symbols are the numerically obtained SD values 
for an uncoupled oscillator. 

is sensitive to the definition of phase. From the observations of circadian periods in mice 
described in the introduction 10|, it is possible that the SCN sends signals to the peripheral 
clocks around the onset of a subjective day. An experimental observation of the checkpoint 
dependence in other biological clocks would be a new source of coupling information. 

We thank Hiroshi Ito for valuable discussions. This work was supported by JSPS KAK- 
ENHI Grant Number 23-11148. 

Appendix A: Proof that the SD is independent of the checkpoint phase in a single 
phase oscillator 

We introduce two checkpoint phases denoted by a and /3. By defining the intervals 
^'^fc"^" — '^k ~'^k ^"^^ ^'^T^^ = ^fe ~ ^fc-i5 ^he oscillation periods observed at a and /3 can be 
decomposed as At^ = At^~^" + At^"*^ and At^ = At^~^^ + Atfl^^"", respectively. Because 
^{t) is independent, the processes a — > /3 and /3 — )■ a for any k are independent. We thus 
have 

E[Atf-"] = E[Atf-1, (Al) 
E[(Atr")2] = i?[(AtJ-")^], (A2) 
10 



and 



El^rf^ Atl-^-] = E[Atr^]E[Atl^n = ^[Atr'AtM. (A3) 
The average and the mean square period are independent of the checkpoint phase labels; 



I.e., 



and 



E[Att] = ^[Atf-"] + E[Atr^] = E[AC,-] + E[Atr'] = E[Att] = r (A4) 



E[{Attr] = E[iAti-^y] + E[iAtr^Y] + 2£;[Aif-"][Atr1 = E[{^4Y]. (A5) 
Thus, we arrive at 



SD(a) = \ E[{Att - r) 



for any arbitrary checkpoint phases a and (3. 



- r?] = SD(/3) 



(A6) 



Appendix B: Calculation of the correlations 



The correlations of the noise terms, ^x,Y{t,4>{t)) = v DZ{(j){t)){E,i{t) ± ^2(^)), are given 



as 



E[U{s, <P{s))U{s\ <P{s'))] = 2DZ(0(s))Z(0(s'))5(s - s'), 
E[^y{s, (Pismris', (Pis'))] = 2DZ{^{s))Z{4>{s'))6{s - s'), 
E[^x{s,<P{s))^Yis',<P{s'))] = 0. 

Using Eqs. (9), ([BT]), and E[^x,y(i, = 0, we obtain 



(Bl) 
(B2) 
(B3) 



E[X{t)^]=exp[2KFx{(Pm 



E[X{Oy' 



+ 



t rt 



^0 



exp[-K{Fxi(Pis)) + Fx{(l){s'))}]E[^x{s, (Pis))^xis' , (l){s'))]dsds' 



= exp[2nFx{(p{t))] 
E[F(t)2]=exp[2/€Fy(0(t))] 



E[X{0)^] + 2D / Z(0(s))2exp[-2fi:Fx(0(s))]rfs 
Jo 

E[Ym 



(B4) 



t rt 



JO 



exp[-«:{Fy + Fy{<P{s'))]]E[^y{s, 0(s))ey(s', (t>{s'))]dsds' 



exp [2fi;Fy(</)(t))] 



E[Y{Q)']+2D I Z((/.(s))'exp[-2KFy(0(s))]ds 





(B5) 
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and 



E[Xit)Yit)] = exp + Fy(0(t)))]E[X(O)F(O)]. 

Substituting t = tcp + t in Eqs. flB4p and flB6p . we obtain 



(B6) 



E[X(tep + r)2] = exp [2kFx( 



yep 



2vr)] 



E[X{0y] + 2D I ' Z(0(s))2exp[-2ftFx(0(s))]rfs 







exp[2/€Fx(ecp)] 



E[X{Of] + 2D 



£cp+T 



Z((/)(s))2exp[-2KFx((/'(s))]rfs 



D[X(t,p)2]+2Dexp[2«^Fx(^cp)] / exp[-2«:Fx(0(s))]c?s, (B7) 

Jo 



and 



E[X(tcp + r)F(t,p + r)] = exp [niFxiO^p + 27r) + Fy(0ep + 27r))]E[X(0)F(0)] 

= expHE[X(t,p)F(tep)], (B8) 

where we use 0(tcp + r) = ^^p + 27r, + 27r) = Fx (6), Fy{0 + 27r) = Fy{0) + c, and 

Z{6 + 27r) = ^(6'). Inserting Oiitcp) = O^p and 0(tcp) = into the definition Ai(t) = 
6*1 (t) — (j){t), we obtain 



Ai(t 



cpV 



0. 



(B9) 



We then obtain 



E[X(t,p)2] 



-E[X{t,p)Y{t,p)] = D[F(tcp)'] = d{9,pf 



(BIO) 



Using Eqs. flBTD- flBTOj) . the relation E[Y{tcp + r)^] = £;[F(t,p)2], and Eq.(lO), we obtain 
the following expression for the phase diffusion 



a{U,p)- = E[{Ai{t,p + r) - Ai{t,p)y 



- {E[X{t,p + rf] + E[Y{t,p + r)2] + 2E[X{t,p + r)F(t,p + r)]} 

^ |2L)exp[2KFx(^cp)] ^ exp[-2KFx(0(t'))]^^^' + 2(1 - exp [Kc])d{9, 



cp) 



C2d{9, 



cpj 



(Bii: 
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Appendix C: Transformation from phase diffusion to period variability 

Here, we illustrate the relationship between o"(^cp) and SD(^cp)- Figure Hl^a) presents a 
schematic view of the trajectories of {/)(t) and 6i{t). An enlarged view of the region around 
(^cp + T, Gcp + 2vr) is displayed in Fig. lU^b), in which the vertical width between the dotted 
lines represents the standard deviation of the phase distribution of 6i(t). In particular, the 
vertical arrow represents the standard deviation of 6*1 (t^p + t), which is denoted by alOcp). 
Because we assume a low noise intensity, the actual trajectories of Oi{t) (thin lines) are very 
close to that of </'(t). We can thus expect that the trajectories are approximately straight 
and parallel to in this enlarged region. Therefore, the horizontal width between the 
dotted lines at 6i = O^p + 2tt is approximately equal to SD(^cp) (horizontal arrow), and the 
relation cr(6'cp)/SD(6'cp) = 0(6'cp) holds approximately. 




FIG. 4. Illustration of the relationship between cj(0cp) and SD(0cp)- 

Appendix D: Definition of phase in the FitzHugh-Nagumo model 

We define the phase ^ as a function of (V,iy), which are the state variables of the 
FitzHugh-Nagumo model given in Eq. (17) in the main text, as follows (Fig. [5]). We first 
assign values to all points on the limit cycle trajectory generated by Eq. (17) without noise 
such that (f) identically satisfies = 27r/r, where r is the period. The limit cycle trajectory 
is independent of the coupling strengths. We set </> = at V = 0.6 with V > 0. We then 
consider radial lines extending from an arbitrary point inside the limit cycle, which we chose 
as (0.6, 0.05) (filled square) in this case. When a radial line intersects the limit cycle at a 
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point that has a value of 0, the phase 9 of all points on the radial line is defined by 6* = 0. 
These radial lines are different from isochrones that give a standard definition of the phase 



18|, but the isochrones are usually unknown. As shown in Fig. 3, our theory is valid even 



for this practical definition. 




FIG. 5. Illustration of the definition of the phase in the FitzHugh-Nagumo model. The limit cycle 
trajectory is generated by a coupled FitzHugh-Nagumo model without noise, whose parameters 
are given in the main text. The circles are placed at equally spaced intervals of (j). All points on a 
straight line radiating from the origin (filled square) have the same phase. 
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